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Abstract 

In this paper we shall introduce a simple, effective numerical method for finding 
differential operators for scalar and vector- valued functions on surfaces. The 
key idea of our algorithm is to develop an intrinsic and unified way to compute 
directly the partial derivatives of functions defined on triangular meshes which 
are the discretization of regular surfaces under consideration. Most importantly, 
the divergence theorem and conservation laws on triangular meshes are fulfilled. 
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Conservation law. 



1. Introduction 

Numerical methods to compute partial differential operators on regular sur- 
faces have always received great interest over last decades. However, they are 
still not well-understood. For example, Conservation laws for diffusion equa- 
tions are usually unsatisfied. Conservation Law is an important principle in 
physics. Indeed, Conservational laws plays a key role in the study of partial dif- 
ferential equations and have many applications in the linearization, integrability 
and numerical analysis. The solution u{x,t) of diffusion equation 

u t = aA^u (1) 

with a > on a regular surface E preserves the total energy. That is, 

^Ju(x,t)dx = 0. (2) 

To numerically simulate u(x,t), one discretizes the regular surface E to ob- 
tain a triangular surface mesh S of E, and considers the discrete solution u(x, t) 
on S. In this way, a fundamental problem arise. Usually, the solution u(x,t) on 
S will not preserve the total energy. That is, 

u(x,U) (3) 
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will change as U increases. 

The violation of the Conservation Law comes from the discretization of the 
Laplacian-Beltrami operator As on E. In this paper, we shall try to handle this 



defect. Lai et. al. [14[ discussed this problem for regular curves in 2008. 

Partial differential equations (PDEs) need to be solved intrinsically and nu- 
merically for data defined on 3D regular surfaces in many applications. For in- 
stance, such examples exist in fluid dynamic flows (Diewald, Preufer and Rumpf 
0), (Bertalmio, Cheng, Osher and Sapiro Q), texture synthesis (Turk 1 18], Witkin 
and Kass[l9|), vector field visualization (Diewald, Preufer and Rumpf [J), weath- 
ering (Dorsey and Hanrahan[10]) and cell-biology (Ayton, McWhirter, McMurty 
and Voth[l[). Usually, regular surfaces are presented by triangular or polygo- 
nal forms. Partial differential equations are then solved on these triangular or 
polygonal meshes with data defined on them. The use of triangular or polygonal 
meshes is very popular in all areas dealing with 3D models. However, it has 
not yet been a widely accepted method to compute differential characteristics 
such as principal directions, curvatures and Laplacians (Chen and Wu|l,Q, Wu, 
Chen and Chi[20|, Taubin[16|). In Wu, Chen and Chi[20[, the authors proposed 
a new intrinsic simple algorithm, LTL method, to handle this difficulty. In this 
note, we shall use this new technique to approximate gradient, divergence and 
Laplace-Beltrami operators on surfaces. 

In Osher and SethianflBt and Bertalmio, Cheng, Osher and Sapiro [2] dis- 
cussed a framework, the implicit surface algorithm, to solve Variational prob- 
lems and PDE's for scalar and vector- valued data defined on regular surfaces. 
Their key idea is to use, instead of a triangular or polygonal representation, 
an implicit representation. The surface under consideration is the zero-level 
set of a higher dimensional embedding function. Then they smoothly extend 
the original data on the surface to the 3D domain, adapt the PDE's accord- 
ingly, and implement all the numerical computations on the fixed Cartesian 
grid corresponding to the embedding function. The advantage of their method 
is the use of the Cartesian grid instead of a triangular mesh for the numerical 
implementation. 

The discretizations of the gradient, divergence and Laplace-Beltrami opera- 
tors that we will discussin this paper will have the following advantages: 

Intrinsicness: we use the intrinsic geometric LTL method to define these op- 
erators. 

Conservation: our Laplace-Beltrami operator will satisfy conservation laws on 
triangular meshes for diffusion equations. 

Convergence: our gradient, divergence and Laplace-Beltrami operators will 
have the linear convergence rate locally and uniformly. 

Simplicity: the numerical computations are also very easy to implement. 

The rest of this paper is organized as follows. In section 2, we recall the 
gradient, divergence and Laplace-Beltrami operators defined on regular surfaces. 
In section three we propose our new discrete algorithm for these differential 
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operators on triangular meshes. We also discuss the convergence problem and 
conservation laws for these operators. Numerical simulations are presented in 
section 4. 



2. The gradient, divergence and LB operators on regular surfaces 

In order to describe the gradient, divergence and the LB operator on func- 
tions or vector fields in a regular surface E in the 3D Euclidean space R 3 , we 
consider a parameterization x : U —> E at a point p, where U is an open subset 
of the 2D Euclidean space R 2 . We can choose, at each point q of x(U), a unit 
normal vector N(q). The map N : x(U) — >• S 2 is the local Gauss map from an 
open subset of the regular surface E to the unit sphere S 2 in the 3D Euclidean 
space R 3 . The Gauss map TV is different iable. Denote the tangent space of E 
at the point p by TE p = {v G R 3 |v _L N(p)}. The tangent space TE P is a linear 
space spanned by {x u ,x v } where u, v are coordinates for U. The gradient V E g 
of a smooth function g on E can be computed from 

^ 9uG-g v F g v E-g u F 

V ^ = EG - F 2 Xu + EG - F 2 x " (4) 
where E,F, and G are the coefficients of the first fundamental form and 

9u = and 9v = (5 ) 

See do Carmo[8] for the details. 

Let X — Ax u + Bx v be a local vector field on E. The divergence, V E • X, 
of X is defined as a function V E • X : E — >• R given by the trace of the linear 
mapping Y(p) — » Vy( p )X for p G E. A direct computation gives 

The LB operator A E acting on the function g is defined by the integral 
duality 

(A E 0,y>) = -(V E 0,V E ¥>) ( ? ) 

for all smooth function ip on E. That is, A^g = V E • V E #. A direct computation 
yields the following local representation for the LB operator A E g on a smooth 
function g: 



VEG-F 2 



y/EG — F 2 



d ( G dg} _ d ( G dg} 
du V y^EG-F 2 du) <9u V ^JEG—F 2 dv ) 

d / G dg\ _ ( G dg_- 
dv v yjEG—F 2 dv' <9tA ^EG-F 2 du > 



(8) 
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3. Discrete gradient, divergence and LB operators 



In this section we shall describe a simple and effective method to define the 
discrete gradient, divergence and LB operator on functions or vector fields on a 
triangular mesh. The primary ideas were developed in Chen, Chi and Wupj lioj 
where we try to estimate the discrete partial derivatives of functions on 2D 
scattered data points. Indeed, the method that we shall use to develop our 
algorithm is divided into two main steps: first we lift the 1-neighborhood points 
to the tangent space and obtain a local tangential polygon. Second, we use 
some geometric idea to lift functions or vectors to the tangent space. We call 
this a local tangential lifting (LTL) method. Then we present a new algorithm 
to compute their gradients in the 2D tangent space. This means that the LTL 
process allows use to reduce the 2D curved surface problem to the 2D Euclidean 
problem. 

Consider a triangular surface mesh S = (V,F), where V = {vi\l < i < ny\\ 
is the list vertices and F = {Tk\l < k < np} is the list of triangles. 

3.1. The local tangential lifting (LTL) method 

To describe the local tangential lifting (LTL) method, we introduce the local 
tangential polygon at a vertex v of V as follows: 

1. The normal vector Na(v) at the vertex v in S is given by 

uotNt 

TeT(v) 

where Nt is the unit normal to a triangle face T and the centroid weight 
is given in [3] by 



_ \\G T -v\\ 2 nn x 
UJ T - — ^ i ( iU J 

„ 2^ \\G t -v\\ 
TeT(v) 

Here, Gt is the centroid of the triangle face T determined by 

Gr = V -±^. (11) 

Note that the letter A in the notation Na(v) stands for the word "Ap- 
proximation" . 

2. The approximating tangent plane TSa(v) of S at v is now determined by 
TS A (v) = {weR 3 \w± N A (v)}. 

3. The local tangential polygon Pa(v) °f v i n TSa(v) is formed by the vertices 
Vi which is the lifting vertex of vi adjacent to v in V. 

Vi = (Vi ~ v)~ <Vi-v, N A (v) > N A (v) (12) 

as in figure [H 
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Figure 1: The local tangential polygon Pa(v) 



4. We can choose an orthonormal basis ei,e2 for the tangent plane TSa(v) 
of S at v and obtain an orthonormal coordinates (x,y) for vectors w G 
TSa(v) by w = xei + ye2- We set ^ = x\e\ + 2/^2 with respect to the 
orthonormal basis ei, 

Next we explain how to lift locally a function defined on V to the local 
tangential polygon Pa(v)- Consider a function h on V. We will lift locally the 
function h to a function of two variables, denoted by h, on the vertices V{ in 
Pa(v) by simply setting 

h(xi,yi) = h(vi). (13) 

and h(0) where is the origin of TSa(v). Then one can extend the function 
h to a piecewise linear function on the whole polygon Pa( v ) m a natural and 
obvious way. 

3.2. A new discrete gradient algorithm 

In this subsection we present a new discrete 2D algorithm for the gradients 
of functions on the 2D domains in the x — y plane and also on triangular surface 
meshes. Given a C 3 function / on a domain Q in the x — y plane with the origin 
(0, 0) G Taylor's expansion for two variables x and y gives 

f(x,y) =/(0,0) + x/ aj (0,0)+ J // y (0,0) 

+ ^/„(0,0) + x^(0,0) + ^/^(0,0) + O(r 3 ) [ ) 

when r = ^ x 2 + y 2 is small. 

Consider a family of neighboring points {xi,yi) G O, i = 1, 2, • • • , n, of the 
origin (0, 0). Take some constants c^, z = 1, 2, • • • , n with J^ILi a f = 1- Then 
one has 

E?=i »0-/(o.o)) 

= (Er=i «^i)/x(o,o) + (Er=i «ii/i)/v(o,o) + ±(£Li «^)/xx(o,o) (15) 
+(Er=i «^iyi)/x»(o, 0) + i(ELi «ii/?)/ ra (o, 0) + o(r 3 ). 

We choose the constants c^, i = 1, 2, • • • , n so that they satisfy the following 
equations: 
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0) Eill a i x iVi = 

One can rewrite these equations in a matrix form and obtain 



x 2 i/2, , x n y n 

9 9 9 

Of* r f* •••••• r Y* 

2/l 5 2/2 5 5 2/n 















) 

















(16) 



Therefore, we have, for these solutions c^, 

n n n 

5>i(/(xi, 2/0-/(0,0)) = (^a i x i )/ K (0,0) + (^a i y i )/ !/ (0 ) 0) + O(r 3 ) (17) 



Choose another solutions i = 1, 2, • • • , n, with ElLi ft! — 1 f° r the linear 
system ([15]) . We have 

(Er=i «^)/,(o, o) + (Er=i «ii/i)/»(o, o) = zr=i v<) - m o)) + o(r 3 ) 

(Er=i ft*i)/*(0, 0) + (ELi ftl/i)/»(0, 0) = Er=i WixuVi) - /(0, 0)) + 0(r 3 ) 

(18) 

If the valence n of the origin (0, 0) is at least 5, we can choose the solutions oti 
and ^ so that the following coefficient matrix is invert ible. 

VEILiA^, E?=ift2/< 

Under these circumstances, we can find the gradient V/(0, 0) = (f x (0, 0), «/y(0, 0)) 
by the relation 

/^(o, o)\ _ f EIU Er=i _1 (T,7=i miffa, 2/0 - /(o, °) A , n(V 2 

- veiLi/^, eiu/W vEr=ift(/(^.^)-mo))y 1 

(19) 

Next we discuss how to approximate the gradient of a function on regular 
surfaces. Let E be a regular surface and S = (V, F) a triangular surface mesh of 
E with mesh size r > 0. Consider a vertex v G V. The local tangential polygon 
Pa(v) of v in Ta(v) is formed by the vertices vi which is the lifting vertex of Vi 
adjacent to v in V. Note that 



(vi - v)- <Vi-v, N A (v) > N A (v). 



(20) 



Choose and fix an orthonormal basis 61,62 for the tangent plane TSa(v) of S 
at v and obtain an orthonormal coordinates (x,y) for vectors w G TSa(v) by 
u> = xei + ye2- We set ^ = xie\ + 2/^2 with respect to the orthonormal basis 
ei,e2- Consider a function ft, on V. We will lift locally the function ft to a 
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function of two variables , denoted by /, on the vertices V{ in Pa(v) by simply 
setting 

f(xi,yi) = h(vi) (21) 

and /(0,0) = h(v) where (0,0) is the origin of TSa(v). In this way, we can 
define the approximating gradient by 

v7 Ah (,A - fEi=i a i x *> T^=i^iyi\~ 1 (YTi=i^i{f{xuVi) -/(o,o))\ . . 
VAhW "lE2=ite. EILi Pan) VEr=ift(/(^>2/i)-/(o,o))J (22) 

where ctj,/^ can be computed from Equations (fT6j) . 

Since the approximating normal vector satisfies N^(v) — Na(v) + 0(r 2 ), one 
can tell from Equations ([T6]) -([T9 ]) and obtain easily the following convergence 
theorem. 

Theorem 1. (^Convergence Theorem 1 ) 

Under the notations as above, one has 

Vvh(v) = V A h(v) + 0(r 2 ). (23) 

3.3. A new discrete divergence algorithm on triangular meshes 

In this subsection, we shall use the divergence theorem to give a discrete ap- 
proximation of the divergence of a vector field X defined on a triangular surface 
mesh S = (V, F). Consider a vertex v G V and let Vj, j = 0, 1, • • • , n be the 
neighboring vertices of v with vo = v n . These vertices Vj&re labeled counter- 
clockwise about the normal vector Na(v). Let Tj be the triangle with vertices 
v,Vj and Vj+\. We define the outer normal vectors n(Tj,Vj) and n(Tj,Vj+i) of 
the triangle Tj at the vertex Vj and Vj+i respectively as follows. See figure [2] 

Consider the lifting vectors w\ and W2 of the vectors Vj+\ — Vj and v — Vj to 
the approximating tangent space TSa(v) of S at the vertex vf 



wi = (vj+i - vj)- < v j+1 - Vj,N A (vj) > N A (vj) 
w 2 = (v - Vj)- < v - Vj,N A (vj) > N A (vj) 



(24) 



Now the outer normal vector n(Tj, Vj) of the triangle Tj at the vertex Vj can be 
defined as 

n(T,,) = <^^1>^-Ihill 2 ^ (25) 
3 JJ || <w 2 ,w 1 >w 1 - |K|| 2 w 2 || 

Similarly, we consider the lifting vectors ws and of the vectors Vj — Vj+i 
and v — Vj+i to the approximating tangent space TSa(vj+i) of S at the vertex 
v j+1 : 

( w 3 = (vj - v j+ i)- < Vj - v j+1 ,N A (v j+1 ) > N A (v j+1 ) , 26 . 
\ w 4 = (v - Vj+i)- < v - v j+1 , N A (vj+i) > N A (vj+i) { ^ 

Note that the outer normal vectors n(Tj,Vj+i) and n(Tj + i, Uj+i) are different. 
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(b) The outer normal vectors n(Tj , v ■ ) and n(Tj , v j+l ) 
Figure 2: The outer normal vectors n{Tj,Vj) and n(Tj,Vj+i) 



Under these notations, we can now define the discrete divergence DivaX of 
a vector field on the triangular surface mesh S by 



Div A X(v) = 



Efc=0 l Tfc l 



E;r 1 fe±F^(2<X(^,n(T i ,^)> 
+2 < X(v j+1 ),n(T j ,v j+1 ) > + < X(v j ),n(T j ,v j+1 ) > ( 2T ) 
+ <X(^ +1 ),n(T,,^) >)] 



where \Tj\ denotes the area of the triangle Tj. 

We can extend the divergence D'ivaX to the whole mesh S piecewise linearly 
in a natural way. Therefore, we have the following lemma. 

Lemma 1. Let X be a vector field on the triangular mesh S. The integration 
of the divergence D'ivaX over S is 



[bw a x= £ 



±Dto A X{vi) \ T i\ 

j£N(i) 



(28) 
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Figure 3: The geodesic triangle T (the blue surface bounded by red curves) n the surface S 
(blue surface) 

Lemma [T] along with the definition of D'yvaX gives the following result. 
Theorem 2. (Discrete Conservation Law 1 ) 

Let S = ( V, F) be a triangular mesh without boundary and X a vector field on 
S. We have 

[ Div A X = 0. (29) 

Let E be a regular surface and S = (V,F) a triangular surface mesh of E 
with mesh size r > 0. If the mesh size r is sufficiently small, we can find a unique 
geodesic 7^ joining two adjacent vertices V{ and vj. In this way, every triangle 
T G F has a corresponding geodesic triangle T on the surface E with the same 
vertices as T. See figure 03 Since the approximating normal vector satisfies 
Nx(v) = Na(v) + 0(r 2 ), the outer normal vectors n(Tj,Vj) and n{Tj,Vj+\) of 
the geodesic triangle Tj in E at the vertices Vj and vj+i respectively also have 
the relations 

f n(T 3l v 3 ) =n(T J ,v J ) + 0(r 2 ) 

\ n(Tj,v j+1 ) =n(T 3 ,v 3+1 ) + 0(r 2 ) V } 

The main purpose of this section is to prove the following result. 
Theorem 3. ^Convergence Theorem 2 ) 

Let E be a regular surface and S = (V, F) a triangular surface mesh of E with 
mesh size r > 0. Consider a smooth vector field X on E ; one has, for sufficiently 
small r > 0, and v G V, 

Div E X(v) = Div A X(v) + 0(r) (31) 
9 



According to the Divergence Theorem on regular surfaces, one has 

Div E X = / <X,n> (32) 
w JdW 

where the domain W is the union of the geodesic triangles Tj with vertices 
v, Vj,Vj+i, j = 0, 1, • • • , n and ft the outer normal vector of W. We denote and 
parametrize the geodesic edge Ej from Vj to Vj+i in Tj by Ej(i), t G [0, 1] with 
= L(^). Then (HQ) gives 

Iw Di vsX = / a < X, n > 

= £"=0% <^,n> (33) 
= Ei=oi(^)/o 1 <^(*)^(*)>* 

where L(Ej) is the length of the geodesic edge Ej. We can approximate the 
vectors X(t) and n(t) by 

X(t) = (l-t)X(v j ) + tX(v j+1 ) + 0(r>) 

n(t) =(l-t)n(f j: v j )+tn(f j: v j+1 ) + 0(r 2 ) 1 ; 

These relations follow from the following easy lemma from Calculus. 

Lemma 2. Consider a smooth function or vector field g on [0,a],a > and a 
sufficiently small r > 0. Then one has, for t G [0,1], 

g(tr) = (1 - t) 5 (0) + tg(r) + 0(r 2 ) (35) 

Equations (JSU) and (J33J) imply 

fj(t) = (1 - t)n(T jjVj ) + tn{T j ,v j+1 ) + 0(r 2 ) (36) 

Note also that the length L(Ej) can be approximated by 

L(E j ) = \\v j+1 -v j \\(l + 0(r 2 )) (37) 

Hence one obtains 

ZU L (Ej)Jo <X(t),n(t)>dt 

= E?=o + 0(r 2 )) [2 < XMMTjW) > (38) 

+2 < X(v j+1 ),n(Tj,v j+1 ) > + < > 
+ <X(v j+1 ),n(T j ,v j ) >]. 

On the other hand, we also have 

J^DivsX =\W\(DivxX(v) + 0(r)) 

= E;=ol^l(Div S X(t;) + 0(r)) (39) 
= E"=o \Tj\(l + 0(r 2 ))(Div s X(«) + 0(r)) 
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Therefore we yield 

+2 < x(v j + 1 ,n(T j ,v j + 1 ) > + < X(vj),n(Tj,v j+1 > ( 40 ) 
+ <X(v j+1 ),n(T j ,v j )>] + 0(r) 
= Div A X(v) + O(r) 

and this proves the main theorem (Theorem [3]). 

Using the results in subsection 3.2 and this subsection, we can approximate 
the Laplace-Beltrami operators on regular surface as follows. Consider a smooth 
function /i on a regular surface E and S = (V,F) a triangular surface mesh 
of E with mesh size r > 0. One can use Equations (fT6]) - (fT9j) to define the 
approximating gradient V A h{v) at a vertex v G V. Then Equation ((6]) gives the 
approximating Laplace-Beltrima operator by 

A A h(v)=Div A (V A h)(v) (41) 

Theorems [T] and [3] then give 

Theorem 4. ( Convergence Theorem 3 ) 

Let E be a regular surface and S = (V, F) a triangular surface mesh of E with 
mesh size r > 0. Consider a smooth function h on S ; one has, for sufficiently 
small r > 0, and v G V , 

Avh(v) = A A h(v) + O(r) (42) 

We can extend A^ft to the whole triangular surface mesh S piecewise linearly 
in a natural way. Then, the Discrete Conservation Law (Theorem [2]) also holds 
for the Laplace-Beltrami operators. 

Theorem 5. (Discrete Conservation Law 2 ) 

Let S be a triangular surface mesh without boundary and h a function on S. 
We have 

[ A A h = (43) 
Js 

Remark 1. The error terms 0{r), 0(r 2 ) in Theorems U\ -0 can be shown to 
depend only on curvatures, injectivity radius fjj [73/ ofH, vector fields X and/or 
functions h. 

Remark 2. We also would like to point out that the methods discussed in this 
section also work in higher dimensions. Namely, we can also use these ideas 
to approximate the gradient, divergence and the Laplace-Beltrami operators for 
hypersurfaces in nD Euclidean spaces with n > 3. We will discuss these in 
another paper. 



4. Numerical simulations 

The Laplace-Beltrami operator on a regular surfaces plays an important role 
on PDEs. In this section, we shall estimate the Laplace-Beltrami operators on 
triangular meshes by our proposed method and shows some numerical simula- 
tions about several important PDEs on regular surfaces. 
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Figure 4: The Laplacian of random polynomial functions on a unit sphere 
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Figure 5: The Laplacian of random polynomial functions on a torus 

4-1. Comparisons of Laplacian estimations 

We compare our proposed method, the level set method and some other dis- 
crization methods for estimating the Laplacian of random polynomial functions 
of degree less than 5 on a unit sphere and a torus in figures 3] and [5] Xu's 
method is a discretization method proposed in 2004. One can find the details 
about Xu's method and Level-set method in [15|, |2lj. We choose 10,000 random 
polynomial functions on these surfaces. The and I2 errors are used for all 
vertices on the triangular mesh. From our simulations, all of these methods are 
convergent and comparable. 

4-2. PDEs on surfaces 

In this subsection, we show numerical solutions of some PDEs on surfaces 
via our proposed method. First, we consider the diffusion equation on a sphere 

u t = A^u (44) 

with the initial condition 

uo(0,rj) = cos(t?) (45) 
where (r, 77) is the spherical coordinate of the unit sphere. We calculate 

u(<j>, 0, t) = exp(-2t) cos(t?) (46) 

as the exact solution of equation ([44]) with initial condition ([45]) . We compute 
the Laplace-Beltrami operator in equation ([44]) by our proposed method 
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Figure 6: The diffusion equation on a unit sphere 



and compare our numerical solution of equation ([44]) with the exact solution 
(pRjj) . Figure [6] illustrates the numerical solutions of equation (jHJ) at time t=0, 
0.5, 10 and 9,000. The "fvals in [a, b] v in the figure [6] means the values u(t) 
on the surface between a and 6, and the "1-infty error" means the error 
of our simulations. Figure 14.21 gives the error of our numerical solutions. 
Obviously, our numerical solution approaches the exact solution when the time 
is large enough. Furthermore, the integration of u, u(t]p)dS, is preserved at 
all time. 

Next, we solve the fourth order diffusion equation, 

ut = -A^A^u (47) 
on the sphere with the initial condition 

u (r, 0, rf) = sin(30)sin(7r]). (48) 



One can find the details about this equation in Greer's paper [12|. In our ex- 
ample, the number of triangles on a triangular mesh is 4096. Figure displays 
the solution at t = 0, t = 0.01, t = 0.5 and t = 5. Obviously, our solution and 
Greer's numerical solution[12] are comparable. 

For our final example, we compute the Allen- Cahn equation, 

u t = e 2 A^u + u 3 -u, (49) 



13 



L 

error 



40 



60 



80 



100 



120 



time t 



Figure 7: The loo error of the diffusion equation on a unit sphere 
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Figure 8: Linear fourth order diffusion on a unit sphere. 
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Figure 9: Allen equation on a torus 



with the initial condition 

f(o,v) 

on a torus, 



-1 otherwise 



(50) 



x(#, rj) = ^(i cos Tj + 1) cos#, cost? + 1) sin 0, ^ sin 77^ . (51) 

Figure [9] shows the results. Again, our results and Greer's numerical solutions 
are equall, well. 



5. Conclusion 

Our proposed method is a new discretization method for estimating the di- 
vergence of a vector field on surfaces. The convergence ratio of our proposed 
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method is as good as the other well-known convergence methods for estimating 
the Laplace-Beltrami operators. Almost all of other methods does not obey the 
divergence theorem, however our proposed method does. That is, our proposed 
method for estimating the Laplace-Beltrami operator on the heat equation have 
the conservation property. In the near future, we shall use our proposed method 
to improve more partial differential equations, such as the Navier- Stokes equa- 
tion, on regular surfaces, triangular meshes and general manifolds of dimension 
n > 3. 
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